/*********************************************************
  SmartMain.cpp
  SMART -- Symbolic Moment Analysis for Rapid Timing
  06-08-2008
  (c) Guoyong Shi;  SJTU
  09-11-2008
  Zhigang Hao; SJTU
 **********************************************************/
//#include  "complex.h"
#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include <cstdlib>
//#include <iomanip>
//#include  <stack>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//#include "globals.h"
#include "utility.h"
#include "ToGraph.h"
#include "RDD.h"
#include "MOR.h"
#include "SmartGlobal.h"
#include "RDDStruct.h"

using namespace std;
// int SampleNum;
// <toGraph> has to be global because it is used by the parser.
class ToGraph *toGraph = new ToGraph;

ofstream MomentFile;
ofstream IndexFile;
FILE  *DelayFile;
FILE  *MOR_File;
FILE  *ORIG_File;
FILE  *Ceff_File;


void Compute_Three_Moment_Delay(double *m1, double *m2, double *m3);
double Search_Delay(double gamma1, double tau1, double tau2, 
    double m1, int nodeN);
void Compute_Effective_Cap(double *m1, double *m2, double *m3);

int main(int argc, char **argv)
{
  char  *ckt_file = NULL;
  char  *moment_file = NULL;
  char  *delay_file = NULL;
  char  *mor_file = NULL;
  char  *orig_file = NULL;
  char  *index_file = NULL;
  char  *ceff_file = NULL;


  int  len;
  int SampleNum;
  char   *stripped_filename; // filename without extension
  const char *moment_suffix = "_moment.out";
  const char *delay_suffix = "_delay.out";
  const char *mor_suffix = "_mor.out";
  const char *orig_suffix = "_orig.out";
  const char *index_suffix = "_index.out";
  const char *ceff_suffix = "_ceff.out";


  class RDD  *rdd; // Resistor Decision Diagram
  class MOR  *mor; // Model Order Reduction
#ifdef HSPICE_SAMPLE
  //ifstream HspiceSample("rlc_tree1000_hspice_mc_var10_sample1000.out",ios::in);
  ifstream HspiceSample("rc_tree5000_hspice_mc_var10_sample1000.out",ios::in);
  if (!HspiceSample)
  {
    cerr<<"file could not be opened."<<endl;
    exit(0);
  }
#endif


  if (argc != 2)
  {
    cout << "\nUsage: SMART netlist.cir \n";
    exit(1);
  }
  // Get the input file name and make an output file name.
  ckt_file = argv[1];

  string s(ckt_file);
  int  dot_at = s.find(".sp");  // the dot position

  // If no extension is found, append ".sp".
  if ( dot_at == -1)  
  {
    stripped_filename = CopyStr(s.data());
    s.append(".sp");
    ckt_file = CopyStr(s.data());
  } else {  // file extension exists
    len = strlen(ckt_file);
    stripped_filename = (char *) malloc(len*sizeof(char));
    strncpy(stripped_filename, ckt_file, len-3);
    //printf("\nstripped_filename = %s", stripped_filename);
  }
#ifdef PRINT_MOMENT
  // (1) Make file name "netlistname_moment.out".
  len = strlen(stripped_filename) + strlen(moment_suffix) + 1;
  moment_file = (char *) malloc(len*sizeof(char));
  sprintf(moment_file, "%s%s", stripped_filename, moment_suffix);
  //s.replace(s.find(".sp"), 11, "_moment.out");
#endif

#ifdef PRINT_DELAY
  // (2) Make file name "netlistname_delay.out".
  len = strlen(stripped_filename) + strlen(delay_suffix) + 1;
  delay_file = (char *) malloc(len*sizeof(char));
  sprintf(delay_file, "%s%s", stripped_filename, delay_suffix);
#endif
#ifdef MOR_COMPUTE
  // (3) Make file name "netlistname_mor.out".
  len = strlen(stripped_filename) + strlen(mor_suffix) + 1;
  mor_file = (char *) malloc(len*sizeof(char));
  sprintf(mor_file, "%s%s", stripped_filename, mor_suffix);
#endif
#ifdef ORIG_MATRIX
  // (4) Make file name "netlistname_orig.out".
  len = strlen(stripped_filename) + strlen(orig_suffix) + 1;
  orig_file = (char *) malloc(len*sizeof(char));
  sprintf(orig_file, "%s%s", stripped_filename, orig_suffix);
#endif
  // (5) Make file name "netlistname_index.out".
  len = strlen(stripped_filename) + strlen(index_suffix) + 1;
  index_file = (char *) malloc(len*sizeof(char));
  sprintf(index_file, "%s%s", stripped_filename, index_suffix);
#ifdef EFFECTIVE_CAP
  // (6) Make file name "netlistname_ceff.out".
  len = strlen(stripped_filename) + strlen(ceff_suffix) + 1;
  ceff_file = (char *) malloc(len*sizeof(char));
  sprintf(ceff_file, "%s%s", stripped_filename, ceff_suffix);
#endif


  cout << "\nInput file name: " << ckt_file;
#ifdef PRINT_MOMENT
  cout << "\nMoment file name: " << moment_file;
#endif
#ifdef PRINT_DELAY
  cout << "\nDelay file name: " << delay_file;
#endif
#ifdef MOR_COMPUTE
  cout << "\nMOR file name: " << mor_file;
#endif
#ifdef ORIG_MATRIX
  cout << "\nORIG file name: " << orig_file;
#endif
  cout << "\nIndexToNode file name: " << index_file;
#ifdef EFFECTIVE_CAP
  cout << "\nCeff file name; " << ceff_file;
#endif
  cout << endl;


  toGraph->Parse_Netlist(ckt_file);
  toGraph->Report_Netlist_Statistics();


#ifdef  TRACE 
  toGraph->Print_All_Elements();
  toGraph->Print_All_Nodes();
#endif

#ifdef PRINT_MOMENT
  // I'll use stream director "<<" to write the moment file.
  MomentFile.open(moment_file);
  if (MomentFile.fail())
  {
    cerr << "\nCan't create file: " << moment_file << endl;
    exit(1);
  }
#endif

  // I'll use stream director "<<" to write the moment file.
  IndexFile.open(index_file);
  if (IndexFile.fail())
  {
    cerr << "\nCan't create file: " << index_file << endl;
    exit(1);
  }
#ifdef PRINT_DELAY
  // I'll use fprintf() to write the delay file.
  DelayFile = fopen(delay_file, "w");
  if (!DelayFile)
  {
    printf("\nCan't create file: %s", delay_file);
    exit(1);
  }
#endif
#ifdef MOR_COMPUTE
  // I'll use fprintf() to write the delay file.
  MOR_File = fopen(mor_file, "w");
  if (!MOR_File)
  {
    printf("\nCan't create file: %s", mor_file);
    exit(1);
  }
#endif
#ifdef ORIG_MATRIX
  ORIG_File = fopen(orig_file, "w");
  if (!ORIG_File)
  {
    printf("\nCan't create file: %s", orig_file);
    exit(1);
  }
#endif
#ifdef EFFECTIVE_CAP
  Ceff_File = fopen(ceff_file, "w");
  if (!Ceff_File)
  {
    printf("\nCan't create file: %s", ceff_file);
    exit(1);
  }
#endif
  toGraph->Create_CTree(); 
  //toGraph->Print_All_Nodes();
  //toGraph->Print_CTree();
  rdd = new RDD(toGraph);
  for (SampleNum = 0; SampleNum < MONTECARLO_SAMPLE; SampleNum++) {
    if (MONTECARLO_SAMPLE > 1) {
#ifdef HSPICE_SAMPLE
      toGraph->Get_Hspice_Sample(HspiceSample);
#else
      toGraph->Generate_Sample();
#endif
      cout << "#Sample" << SampleNum+1 << endl;
    }
    for (int i = 1; i <= rdd->Get_MaxOrder(); i++) {
      rdd->Evaluate_Moments(i);
    }

    // rdd->Print_Symbolic_Moments();
#ifdef PRINT_MOMENT
    //rdd->Dump_Moments(2, OutFile );  // up to order 2
    //rdd->Dump_Moments(3, cout);
    rdd->Dump_Moments(rdd->Get_MaxOrder(), MomentFile, SampleNum);
#endif
    if (SampleNum == 0)
      rdd->Dump_IndexToNode(IndexFile);

#ifdef EFFECTIVE_CAP
    Compute_Effective_Cap(
        rdd->Get_Moment(1),
        rdd->Get_Moment(2),
        rdd->Get_Moment(3));
#endif

#ifdef PRINT_DELAY
    Compute_Three_Moment_Delay(
        rdd->Get_Moment(1),
        rdd->Get_Moment(2),
        rdd->Get_Moment(3));
#endif

#ifdef MOR_COMPUTE
    mor = new MOR(rdd->Get_Moment_Length(),
        4, // reduced order 
        rdd);

    if (toGraph->Get_Num_Inds() > 0) {
      mor->Gram_Schmidt_RLC();
      mor->Reduce_C_G_RLC(rdd->Get_CTree());
    } else {
      mor->Gram_Schmidt_RC();
      mor->Reduce_C_G_RC(rdd->Get_CTree());
    }

    //mor->Print_V(stdout);
    mor->Print_V(MOR_File);
    //mor->Print_Reduced_V(MOR_File);

    //mor->Print_Reduced_C_G(stdout);
    mor->Print_Reduced_C_G(MOR_File);
#ifdef ORIG_MATRIX
    mor->Print_Orig_C_G(stdout);
    mor->Print_Orig_C_G(ORIG_File);
#endif
    delete mor;
    //mor = 0;
#endif
  }
#ifdef MOR_COMPUTE
  fclose(MOR_File);
#endif

#ifdef ORIG_MATRIX
  fclose(ORIG_File);
#endif

#ifdef EFFECTIVE_CAP
  fclose(Ceff_File);
#endif

#ifdef PRINT_MOMENT
  MomentFile.close();
#endif
  cout << endl;  
}

void Compute_Effective_Cap(double *moment1, 
    double *moment2, 
    double *moment3)
{
  int MomentLength;
  double Res;
  double Cap1;
  double Cap2;
  double alpha;
  double beta;
  double Ceff;
  double m1,m2,m3;
  MomentLength = toGraph->Get_Num_Caps();
  for ( int i = 0; i < MomentLength; i++) {
    //    if (i==1650 || i==2808 || i==848) {
    m1 = moment1[i];
    m2 = moment2[i];
    m3 = moment3[i];
    Res = m3*m3/(m2*m2*m2);
    Cap1 = m2*m2/m3;
    Cap2 = m1-m2*m2/m3;
    alpha = Cap2/(Cap1+Cap2);
    beta = R_EFF/Res;
    Ceff = (3*alpha+beta*beta)*(Cap1+Cap2)/(3+beta*beta);
    fprintf(Ceff_File, "\n%-8.2e", Ceff);
    //fprintf(Ceff_File, "\n%-12.2e%-12.2e%-12.2e%-12.2e%-12.2e%-12.2e%-12.2e%-12.2e%-12.2e", m1, m2, m3, Res, Cap1, Cap2, alpha, beta, Ceff);

    //printf("\n%-12.2e", Ceff);
    //    }

  }
}
// Use the first three moments to extract delay.
// Delay calculated only for moments satisfying:
//  gamma = m1*m1/m2 < 1;   (near-end condition)
// Under the condition (gamma < 1), two positive time constants
// tau1 & tau2 are determined. 
// Then, tau1, tau2, and m1 are used to find a 50% delay via a binary search.
//
void Compute_Three_Moment_Delay(double *moment1, 
    double *moment2, 
    double *moment3)
{
  // Length of moment vector
  int   LenMom;
  char  *node_name;
  double m1, m2, m3;  // three moments at one node
  double m3_2, m2_1; // m3_2 = m3/m2, m2_1 = m2/m1
  double  gamma1, gamma2;  // gamma2 = gamma1^2
  double A, B, Delta, Delta1;
  int  positiveDelta = 1;
  double tau1, tau2;  // two time constants
  double ED, D2M;  // Elmore Delay and D2M delay
  double delay;   // the delay estimated using tau1, tau2
  double delay_scaled; // the estimated delay scaled by gamma1

#ifdef  CONSOLE
  printf("\n%-8s%-12s%-12s%-12s%-10s%-12s%-12s%-12s",
      "Node", "Elmore", "D2M", "gamma1", "A > 0", "tau1", "tau2", "delay~");
#else
  fprintf(DelayFile, "\n%-8s%-12s%-12s%-12s%-10s%-12s%-12s%-12s",
      "Node", "Elmore", "D2M", "gamma1", "A > 0", "tau1", "tau2", "delay~");
#endif

  // fprintf(DelayFile, "\n%-8s%-12s%-12s%-12s%-12s%-12s",
  //   "Node", "Elmore", "D2M", "gamma1", "delay", "delay_scaled");


  LenMom = toGraph->Get_Num_Caps();

  // Process all nodes in order.
  for ( int i = 0; i < LenMom; i++) {

    m1 = moment1[i];
    m2 = moment2[i];
    m3 = moment3[i];

    m3_2 = m3/m2;
    m2_1 = m2/m1;

    gamma1 = m1 / sqrt(m2);
    gamma2 = gamma1*gamma1;

    // Compute two delay estimates
    ED   = log(2)*m1;   // Elmore Delay
    D2M = log(2)*m1*gamma1; // D2M delay

    // Compute two time-constants using three moments
    A = (m3_2-m2_1)/(m2_1-m1);
    Delta1 = (A+1)*(A+1)-4*A*gamma2;

    if (Delta1 >= 0) {
      B = sqrt(Delta1);
      tau1 = 0.5*m2_1*((A+1)+B);
      tau2 = 0.5*m2_1*((A+1)-B);
    } else {
      printf("\nCompute_Three_Moment_Delay(): Negative Delta!\n");
      positiveDelta = 0;
    }

    node_name = toGraph->Get_Node_By_Index(i)->name;

    // For gamma1 < 1, it is possible that A < 0, which implies
    // that one of the time constants is negative.
    //if (gamma1 < 1 && A > 0) {
    if ( A > 0 && positiveDelta ) {
      delay = Search_Delay(gamma1, tau1, tau2, m1, i);
      //delay_scaled = delay * gamma1;

#ifdef CONSOLE
      printf("\n#%-7s%-12.2e%-12.2e%-12.3f%-10.3f%-12.2e%-12.2e%-12.2e", 
          node_name, ED, D2M, gamma1, A, tau1, tau2, delay);
#endif

#ifdef  PRINT_TAU   
      fprintf(DelayFile, "\n#%-7s%-12.2e%-12.2e%-12.3f%-10.3f%-12.2e%-12.2e%-12.2e", 
          node_name, ED, D2M, gamma1, A, tau1, tau2, delay);
#else
      fprintf(DelayFile, "\n#%-7s%-12.2e%-12.2e%-12.3f%-12.2e%-12.2e", 
          node_name, ED, D2M, gamma1, delay, delay_scaled);
#endif
    } else {
      // Delta is negative; we use D2M
      delay = D2M;

#ifdef CONSOLE
      printf("\n#%-7s%-12.2e%-12.2e%-12.3f%-10.3f%-12s%-12s%-12.2e", 
          node_name, ED, D2M, gamma1, A, "----", "----", delay);
#endif

#ifdef  PRINT_TAU   
      fprintf(DelayFile, "\n#%-7s%-12.2e%-12.2e%-12.3f%-10.3f%-12s%-12s%-12.2e", 
          node_name, ED, D2M, gamma1, A, "----", "----", delay);
#else
      fprintf(DelayFile, "\n#%-7s%-12.2e%-12.2e%-12.3f%-12.2e", 
          node_name, ED, D2M, gamma1, delay);
#endif
    }
#if 0  
    printf("\n#%s\t%.2e\t%.2e\t%.3f\t%.3f\t%.2e\t%.2e\t%.2e", 
        node_name, ED, D2M, gamma1, A, tau1, tau2, delay);
#endif

  }
  fprintf(DelayFile, "\n");
  fclose(DelayFile);
  cout << endl;
  }
  // Search (via bisection) the 50% delay (t_{0.5}) of the step response
  //  y = 1 - lambda1*exp(t/tau1) - lambda2*exp(t/tau2).
  // t_{0.5} is supposed to lie in the interval [0, m1].
  double Search_Delay(double gamma1, double tau1, double tau2, 
      double m1, int nodeN)
  {
    double lambda1, lambda2;
    double T1, T2, t, y;
    int  N;
    int  loopLimit = 10;
    char  *node_name;

    lambda1 = (m1-tau2)/(tau1-tau2);
    lambda2 = (tau1-m1)/(tau1-tau2);

    // Do the first check (gamm1 near to 1)
    t = (tau1+tau2)/2;
    y = 1 - lambda1*exp(-t/tau1) - lambda2*exp(-t/tau2);
    if (fabs(y-0.5) < 1e-3) {
      return (t); // tau1 and tau2 are close to each other
    }
    if (gamma1 <= 1) {
      T1 = 0;
      T2 = m1*log(2); 
      t = T2;
    } else {
      // Find an interval that includes  t_{0.5}
      T1 = m1*log(2);
      T2 = 2*T1;
      y = 1 - lambda1*exp(-T2/tau1) - lambda2*exp(-T2/tau2);
      while (y < 0.5) {
        T1 = T2;
        T2 = 2*T1;
        y = 1 - lambda1*exp(-T2/tau1) - lambda2*exp(-T2/tau2);
      }
    }
    N = 0; // counter
    //printf("\n#%d:  y(%.2e) = %f", N, t, y);
    // Bisection search t_{0.5} in the interval [T1, T2]
    while (fabs(y-0.5) > 1e-3) {
      t = (T1+T2)/2;
      y = 1 - lambda1*exp(-t/tau1) - lambda2*exp(-t/tau2);
      if (y > 0.5)   T2 = t;
      else    T1 = t;
      N = N+1;
      //printf("\niter #%d:  y(%.2e) = %f", N, t, y);
      if (N > loopLimit) {
        node_name = toGraph->Get_Node_By_Index(nodeN)->name;
        printf("\nSearch_Delay() Loop limit exceeded at node [%s].",
            node_name);
        break;
      }
    }
    //printf("\ny(%.2e) = %f\n", t, y);
    return  (t);
  }
